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1. Program Summary 

• Title of the program: Hamevol 

• Version number. 1.0 

• Avaible at: http://wwwteor.mi.infn.it/ antonell/programs/RKutta , and mirrors 

• Programming language: C++ 

• Platform: any platform supporting a C++ compiler (examples: Linux, Unix, Windows) 

• Tested on:Pentium PC, AMD PC 

• Memory requirements for execution: Standard application: 40 Kbytes 

• No. of bytes in distributed program, including test data, etc: 235.000 

• Keywords: Numerical algorithms, Differential equations, Hamiltonian evolution. Oscilla- 
tion, 



• Nature of physical problem: Numerical solution of Hamiltonian differential equations. Ap- 
plication to the numerical calculation of the oscillation probability for a quantum system 
(like, for instance, neutrinos of any kind propagating in a medium) 

• Method of solution: Algorithm based on fifth order semi-implicit Runge-Kutta method 

• Typical running time: ~ 30 seconds for every single point on a Penthium IV PC 



2. Introduction. The mathematical problem. 

The use of numerical algorithms and suitable computational techniques has often been very 
useful to find the solution of difficult problems which are of interest for mathematics and other 
applied science. This is particularly true in our days, when the relationships between information 
technology and other sciences are becoming closer and closer. 

In this paper we discuss the adaptation of well known numerical techniques to a general 
class of problems which are described by ordinary differential equations and we present some 
examples taken from physics. 

The numerical code and algorithm we are presenting in this work is based on the imple- 
mentation of the Runge-Kutta method and it can find significant applications in the study of 
different physical systems. In fact the evolution of every system can be fully described once 
we are able to solve the differential equations that drive this evolution. The typical example 
is the case in which one wants to study a linear quantum system that is described by a vector 
X = {Xi,i = 1..N), where Xi are the elements of an appropriate basis decribing the system. 
The system of linear differential equations we are interested in can be put in the simple form: 

i^ = H{t)X, (2.1) 

where H is a matrix determining the evolution of the system. 

In the language of physics H is the Hamiltonian of the system and Equation (|2.l| ) is the 
corresponding Schroedinger equation. It is clear, however, that the system of differential equa- 



tions (2J) is very general and they can describe also problems of different nature in fields that 
are completely different from physics. Depending on the kind of problems one has to solve, 
the requirements which are fundamental for the solutions can be different. This can suggest 
the choice of a particular algorithm in order to fulfill these requirements in the best possible 
way. For instance, as we are going to discuss later, in the physical problems we are intersted in, 
the delicate point is efficiency more than accuracy and this justifies the choiche of a particular 
version of adaptive Runge-Kutta method that, if properly adapted to our purpose, enable us to 
obtain satisfactory results. 

The solution of equation ( |2.lD can be written in terms of the fundamental system of solutions. 



or equivalently called evolution operator ?7(t, ig); defined by the expression: 

U{t,to)=e^p[-iH{t-to)], (2.2) 

where to is the initial time at which we know the state of the system. The simple formula given 
above is valid for a time independent Hamiltonian. In the case in which the Hamiltonian of the 



system is changing in the time (hke, for instance, in presence of matter effects), formula ( p.2| ) 
must be replaced by path-ordered exponential 



U{t,to) = Pexp 



i [ drHlr] 

J to 



(2.3) 



The code we are presenting here can be used to solve in an iterative way the system of 
differential equations appearing in Equation (|2.1| ). This is particularly important in the case of 
Hamiltonians which are explicitly time dependent or which contain terms fastly oscillating in 
time. 

The structure of this paper is the following. We start discussing in section |3| some concrete 
examples, taken from physics, of situations in which the application of the algorithm presented 
here is particularly suited. In the next section we present the algoritmic structure and salient 
aspects of the code we developed. This section includes a presentation of the algorithm, all the 
essential informations about the distribution, the main subroutines and functions and about the 
way of working of a sample program. In section ^ we draw our final conclusions. An example of 
a sample program, with some specific Hamiltonian set up and the relative outputs are presented 
in the Appendices. 



3. The Physical motivations. 

A very interesting example is given by the study of neutrino physics, which has been one of the 
central topics of elementary particle physics in the last years. A detailed discussion about the 
main properties of neutrinos and the relevance of their study for our knowledge of the intimate 
structure of matter is beyond the scope of this work. Therefore we refer the interested reader 
to the many reviews one can find in literature [Q]. Here we just recall that during the last years 
the answer has been given to the central question of neutrino physics, which puzzled physicists 
for more than seventy years, that is to discriminate whether this particle is massive or massless. 

We know by now that neutrinos are massive and oscillating particles and the proof of this 
has been given by the important results obtained mainly in the last decade (and especially in 
the very last years) by the experiments looking for oscillation signals of neutrinos from different 
sources: solar ||2|, ^, reactor [Q] and atmospheric Q neutrinos. 

The great relevance of these results is confirmed by the fact that this is up to now the strongest 
indication of oscillation we have in the leptonic sector and it is impossible to accomodate it in the 
usual "minimal version" of the Standard Model (the theory describing very well the electroweak 
interactions of elementary particles) . One can say that neutrino oscillations is a hint for physical 
phenomena beyond to what is presently known. 

All these experimental evidences in favors of the oscillation hypothesis have proved that the 
flavor eigenstates of neutrinos, that is the ones entering in the weak processes, are in fact quan- 
tum superspositions of different mass eigenstates (at least three different mass eigenstates are 
needed to explain the full set of experimental data, if one doesn't take into account the contro- 
versial results of the LSND experiment). During the evolution the composition of this quantum 
system can change giving rise to the oscillation phenomenon detected in the experiments. 

Any study of neutrino oscillation is necessarily based on the calculation of the so called 
neutrino survival (or transition) probability. This is the probability that a neutrino emitted by 



a source with a certain flavor, for instance an electronic neutrino emitted in the solar fusion 
processes, remains with the same flavour (or is converted into a different flavor neutrino) before 
reaching the detector. 

Hence the basic quantity to compute is the survival probability in matter for a neutrino of 
a certain flavor: 

P{iyi^Ui;t,E^) = \{iyi{to)\U{t,to)\iyi{to)) P (3.1) 

where to is the inital time at which the neutrino is assumed to be in the flavour eigenstate Vi 
(with i = e, fi, r) and U {t, to) is the evolution operator given by Equation (|2.3|). 

In absence of matter the basis in which the Hamiltonian is diagonal is simply the mass basis, 
whose eigenstates Va are connected to the neutrino flavor eigenstates by means of the mixing 
matrix U : 

Vi = Uia Vol-, i = 6, //, t and Q = 1, 2, 3 . (3.2) 



Using the same compact notation of Equation ( p.lD , we denote the set of the three neutrino 
mass eigenstates with the vector where v = (z^o,a = 1,2,3). This notation can be simply 
extended to the case in which one has more than three neutrinos ^ . The Schroedinger equation 
describing the evolution of the system in vacuum under the relativistic approximation and in 
the hypothesis of equal momentum for different mass eigenstates is: 

= }i% (3.3) 

dt ^ ' 

where 

= DiagEa = \/ + ^Ey^\ m\l2E^ (3.4) 

\ mll2E^ ) 

The last experimental data (both for atmospheric and solar neutrinos) have proved that to 
describe neutrino evolution one has to take into account also the modification of the oscillation 
pattern due to the very important effects of interaction with matter. This gives rise to the well 
known MSW effect Q. The problem of calculating neutrino oscillation probability in presence 
of matter effects has been faced by many authors with different approaches, both numerical and 
analytical, in the case of two neutrino flavors @, |8|, |9|. Exact solutions to the three neutrino 



MSW equations were derived |1C, 11] for simple matter densities. Numerical algorithms for 
direct computation of the solar neutrino survival probability with all three active neutrinos have 
been presented in [12, 13, 14| . 



The purpose of our code Hamevol is to calculate the electron-neutrino survival probability 
for a given neutrino energy, a given density profile, given neutrino masses and mixing matrix. 
This set of (three) oscillation probabilities will be used subsequently by other (Fortran) programs 
to calculate expected signals from diverse neutrino experiments. Due to the fact that both the 
mixing matrix and the density function are considered input parameters of Hamevol , all kinds 
of neutrino oscillation problems may be tackled, including those relevant to anti-neutrinos only. 

In presence of standard matter with arbitrary electron number density, the propagation is 
usually well described by the following system of differential equations: 



^This is the case if one introduces also sterile neutrinos in the analysis 
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(3.5) 



where V is a matrix with Vn as only non-zero element, p{t), essentially a forward scattering 
amplitude, is proportional to the electron number density of the medium Ne{t) 

p{t) = ±V2GFNe{t) (3.6) 

and U is the mixing matrix connecting the neutrino flavor eigenstates with the mass eigenstates. 



In the case of three neutrino generations, adopting the Particle Data Group [|15[ convention for 
the mixing matrix, one gets: 
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(3.7) 



where q = cos^j and Sj = sin^j and the three mixing angles 9i = 9 12, O2 = 6*23, and ^3 = C13 
roughly measure mixing between mass eigenstates (1-2), (2-3), and (1-3) respectively. We have 
neglected the CP violating phase, which is irrelevant in this problem. 

The plus/minus sign in formula of Equation ( |3.6| ) is for neutrinos/antineutrinos respectively. 
The time dependence of the electron density Ne is a crucial factor involved in solving the 
evolution equation. 

The difference of the eigenvalues of ( the inverse of the usually defined oscillation length) 
is typically considered to be, in the solar case, of the order 
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an analytical solution of the problem has been found for the case in which the electron 



number density is parametrized (data taken from [16|) , for sufficiently far distances from the 
solar core as: 



Ne{r) = No exp(-Ar) ; A ~ 10.6 ■ 



(3.8) 



with r, ro the distance from the center and solar radius respectively. 

An important peculiarity of neutrino propagation in matter with respect to vacuum is that 
for a certain value of the parameters a resonance appears at a certain point along the neutrino 
trajectory. This resonance takes place if the following condition is satisfied: 



Am^ cos 29 
E '■ 



(3.9) 



where 9 is the mixing angle between the 2 neutrino flavors taking part to the oscillation phe- 
nomenon. Apart from this approximate results, exact solutions have appeared for particular 
forms of the function p: for linear densities, in terms of Weber-Hermite functions (|17, 18|); for 
functions of the form p{t) = C{1 + tanh(At)) in |l^, and for exponentially decaying densities 
p(t) = ce-^* in0,g. 

The parametrization of Equation ( |3.8D is the one used also in the sample program we present 



as an example in subsection 4.4. 



In any case our numerical alghorithm enables us to find a solution of the problem for every 
expression one chooses for the electron density number. 

The capability of solving the system of Equations ( p.5| ) and, therefore, of computing the 
neutrino survival (or transition) probabilities as a function of the mixing parameters is an 
essential ingredient in every analysis of neutrino data. The theoretical expected signal in every 
experiment is obtained by convoluting neutrino fluxes, oscillation probabilities, neutrino cross 
sections and detector energy response functions. The comparison of these expected signal as a 
function of the mixing parameters with the experimental results is then usually performed by 
means of a statistical analysis. The outcome of this kind of analyses is typically the production 
of exclusion plots selecting the regions of the mixing parameter plan which are in agreement 
with the data at a certain confidence level. The algorithm we are presenting in this work enables 
us to numerically solve the neutrino evolution equations for all the oscillation parameter space, 
without the need to introduce the approximation which are required in different approaches 
based on the use of semi-analytical expressions in portions of the parameter space. For a more 
detailed description of the full procedure we adopted for a phenomenological analysis of solar 
and reactor neutrino data we refer the interested reader to . 

The numerical alghoritm we are presenting here finds also other relevant applications. Some 
significant examples are: 

a) the study of the neutrino evolution inside stochastic media and neutrino propagation in solar 
magnetic field |2^]. In these cases one has to solve systems of differential equations where, 
respectively, (2A^)^ and (2A^) equations appear (N is the number of neutrino species); 

b) solution of the renormalization group equation for the Minimal Supersymmetric Standard 
Model in a Supergravity scenario. Here the number of differential equations in the systems that 
have to be solved simultaneously (with a complicated mixture of initial and boundary conditions) 
is typically between thirty and one hundred. 



4. Code Structure 
4.1 Algorithmics 

There are different ways of solving a possibly stiff set of equations and the advantages and 
drawbacks of each of them must be evaluated keeping in mind the kind of problem one wants to 
study. In our case we aim for efficiency rather than precision (a relative precision of lO-^-lO"^ m 
conversion probabilities for example would be sufficient in a typical application in particle physics 
propagation problems). Therefore we opted for an adaptive Runge-Kutta (RK) algorithm. 



The Runge-Kutta method |23] is particularly suitable for solving differential equations, 
starting from the knowledge of the function at the a fixed intial point X and advancing the 
solution from X to X + /i, by using the evaluation of the function at intermediate points inside 
the interval h . By properly combining these evaluations one can reduce the error in the final 
output. The method is conventionally denoted of order n if its error term is O (/i""^^) 

An essential characteristic of a good Ordinary Differential Equations integrator is the capa- 
bility of having an adaptive control over its progress and a mechanism for adapting its stepsize, 
in such a way to obtain the required accuracy with the minimal possible computational effort. 
This property is possessed by the algorithm we implemented. Although an implicit RK would 
be advised for stiff equations, there are several alternatives, such as the semi-implicit fifth-order 



RK routine we have chosen. This routine requires the determination of the function at five 
different points in the interval between the chosen steps. They are 



ki = hf {Xn,yn) 

k2 = hf {xn + a-ih, yn + b2iki) 

(4.1) 

ke = hf {xn + aeh, y„ + ftei^i + ■ • • ^es^s) 

6 

Vn+l = yn + Cjfcj 
1=1 

where the coefficients Oj, bij and q must satisfy certain constraints in order to ensure stabihty 
and convergence. This algorithm is well suited for adaptive stepping, due to the fact that among 



the six evaluations in Eq. (4.1) there is an embedded fourth-order combination, which, although 
redundant, gives us an estimate of the error at each evaluation thereby allowing us to adjust the 
step size. Table (|l|) gives a list of Oj, bij and q as determined by Cash and Karp [23|. 
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Table 1: Cash-Karp coefficients for our RK routines, taken from Rcf.[E3[ 



4.2 The distribution 

The distribution of the program is contained in the tarred gzipped file hamevoll.O.tar.gz. In 
any Linux or Unix system, unpacking and untarring this distribution file will produce a local 
directory called Hamevoll.O containing the following ascii files: 



hamevol-rungekutta.hpp The file with the main routines. 



hamevol-util .hpp 
hamevol-sample . lipp 
hamevol-sample . cpp 
Makefile 



Some auxiliary utilities. 

The sample program header. 

The sample program. 

A simple compiling make example. 



In the same Linux or Unix system the execution of the command 
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. /make 



or directly the explicit call to the C++ compiler 

g++ -03 -o hamevol.x hamevol-sample . cpp 

should produce an executable file from our sample program dedicated to solve some par- 
ticular neutrino propagation problem. To run the executable the following command should be 
typed 

./hamevol.x OPTION 

where in our sample program OPTION=1,0 depending on whether full Sun+Earth or only 
Sun propagation is demanded. As a consequence, some brief ouput will appear on the standard 
ouput, mainly information about parameter settings and options. In addition, as a result of 
the execution our sample program writes an output file runge . out with the neutrino conversion 
probabilities along the neutrino trajectory. 

In our distribution, we have well separated the code corresponding to the general routines 
implementing the Runge-Kutta algorithm from those corresponding to a particular application 
(the "sample" files) of interest to us and that are presented here: the computation of oscillation 
neutrino probabilities. Within our sample programs, it is also well differentiated the driver code 
which calls the RK routines from the part where a concrete hamiltonian is built. In the most 
basic case, a general user should be able to use our program as a black box for his own purposes 
simply plugging his own definition for the hamiltonian. 

In the following we will first describe the main routines included in file hamevol-rimgekutta.hpp 
and then those contained in the sample programs. 

4.3 The RK algorithm. Main Subroutines 

The kern code is built up with five main subroutines. They correspond to procedures and meth- 
ods well known in the literature. We have improved and adapted them for our purposes. The 
following classes are located in file hamevol-rungekutta.hpp. Here it follows a brief description of 
any of them together with its calling sequence. For brevity, arguments which are also referenced 
somewhere else are omitted here. 

• void runge (CNumber y[], CNumber dydx [] , CNumber (*H) (...), int n. Number x. 
Number h., CNumber yout [] , void (*derivs) ( . . . ) ) 

Given the value y [1 , . . ,n] of the vector state describing a phyisical system made up with n 
components and evolving according to Schroedinger equation and knowing the Hamiltonian 
H of the system, the subroutine produces the advanced solution as the function at the 
incremented variables yout [1 . .n] . 

• void ode int (CNumber ystart[], int nvar. Number xl. Number x2. Number eps, 

Number hi, Number hmin, CNumber (*H) (...), void (*derivs) ( . . . ) , void (*rkqs) ( . . . ) ) 
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This is a Runge-Kutta driver with adaptive stepsize control. It Integrates the starting 
values ystart[l. .nvar] from xl to x2 with accuracy eps, storing the intermediate results 
in global variables. 

A value hi should be set as a guessed first stepsize, hmin as the minimum allowed stepsize 
(it can be zero) . On output nok and nbad are the number of good and bad (but retried and 
fixed) steps taken, and ystart is replaced by values at the end of the integration interval. 

• void rkqs (CNumber y[], CNumber dydx [] , int n, Number *x, Number htry, Number eps, 
CNumber yscal[], Number *hdid, Number *hnext , CNumber (*H) (...), void (*derivs) (...)) . 

This routine rkqs is implemented in order to perform an adaptive 5th order Runge-Kutta 
integration. The method enables to have a monitoring of local truncation error, in order 
to ensure the required accuracy and adjust the stepsize. The inputs are the independent 
variable vector y[l. .n] and its derivative dydx[l. .n] at the starting value of the indepen- 
dent variable x. Other inputs are the stepsize htry, the required accuracy eps, and the 
vector yscal[l. .n] against which the error is scaled. On output, y and x are replaced 
by their new values, hdid is the stepsize that was actually accomplished, and hnext is the 
estimated next stepsize. 

• void rkck (CNumber y[] , CNumber dydx[] , int n. Number x. Number h, CNumber yout[] , 
CNumber yerr [] , CNumber (*H) (...), void (*derivs) (...)) 

Used in adaptive size Runge-Kutta integration. Given values for the variables y[l. .n] 
and their derivatives dydx[l. .n] known at x, advance solution over an interval h and 
return the incremented variables as yout[l. .n]. Also returns an estimate of the local 
truncation error in yout using embedded fourth-order method. 

The user supplies the routine H(x,i,j), which returns the element (i,j) of the hamilto- 
nian of the evolution. 

• void deriv(Number x, CNumber y[], CNumber dy[], int n, CNumber (*H) (...)). 

The user-supplied routine derivs is used for calculating the right-hand side derivative. 
The user supplies the routine derivs (x, y, dydx, n,H), which returns the derivatives dydx of 
the many variable function y with respect to the vector x at the point x. 

The following Auxiliary functions which will allocate the following data structures are in- 
cluded in the auxiliary distribution file hamevol-util .hpp. 



template <class Type> Type **matrix( . . . ) allocate a Type matrix with a subscript range. 



CNumber *Cvector (long nh): 
Number *vector(long nh): 
int *ivector (long nh) : 

unsigned long *lvector (long nh): 
unsigned char *cvector (long nh): 



allocate a CNumber vector with subscript range v[l. .nh] 



a vector with subscript range v[l. .nh] . 
a vector with subscript range v[l. .nh] . 
a vector with subscript range v[l . .nh] . 
a vector with subscript range v[l . .nh] . 
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4.4 Sample Program and Inputs 

The void main routine in hamevol-sample . cpp file takes the values of the neutrino wave functions 
at initial starting points for physical initial conditions which are standard for solar neutrino 
physics (fe(0) = 1, f^(0) = 0, ^'r(O) = 0) and calculates the final wave function and corresponding 
probabilities at the target final points. The algorithm includes the following steps: 

• It takes from the command line the user argument "1" (Sun) or "0" (Earth) propagation. 

• It performs argument validation, set internal flags and writes to the standard output a list 
of current values of diverse parameters. 

• It declares the output file stream out ''runge.dat", class of stream included within 
<f streciin.h>. 

• It declares the Cvector objects nu.dnu, respectively instances of the neutrino wave function 
and their vector derivative. It performs diverse other initializations. 

• Finally functions odeint and evolve are repeteadly called until the desired final point or 
the maximum number of steps is reached. 

Different parameters, for example the number of equations, or in this physical case the 
number of neutrino species (A^ = 2,3), have to be set in the header file hamevol-sample .hpp . 
We run the RK algorithm to obtain the transition probabilites of neutrinos produced at the sun 
center with a given energy and mixing angle as a function to its position along the trajectory 
sun-earth. The user should provide routines for computing the electron densities at Sun and the 
Earth along the neutrino trajectory. They appear in the matter part of the neutrino hamiltonian. 
We include examples of the main program and other smaller routines as appendices. 

5. Conclusions 

A new code based on a semi-implicit fifth order adaptive Runge-Kutta algorithm has been 
developed by us. It can be used as solver for many systems of differential equations, like, for 
instance, the ones that usually describe the evolution of a system in physics and in other fields. 
This algorith is particularly suited for the solution of differential equations in which the operator 
driving the evolution of the system is changing in time. 

Here we focus our attention to the application of this code to the study of physical problems, 
like solving the Schroedinger equation for a system that is a quantum superposition of different 
possible states. The explicit example we present is the study of the evolution and calculation 
of transition probabilty for neutrinos emitted by a source and travelling in a medium. This 
code has been already applied by us as a useful tool to obtain a check with respect to other 
possible numerical algorithms (like the ones based on the evolution operator formalism) in 
our phenomenological analysis of different neutrino oscillation experiments. This analysis has 
confirmed the validity of neutrino oscillation hypothesys and enabled us to determine the allowed 
region for mixing parameters, a topic of great relevance in Elementary Particle Physics. 

In this paper we discuss the structure of the algorithm we developed and the main features 
of our code. We also present a sample program and give some typical outputs, as a concrete 
example of application of our algorithm. 
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7. Appendices 
7.1 Using Hamevol 

Here we present the main sample program, corresponding to tlie file hamevol-sample . cpp, wiiere 
we siiow explicitly the use of the main routines. 

#include "hamevol .hpp" 
#include <string.h> 

/* Inizialize the vacuum values ( )*/ 

/* #define parameters ( )*/ 

Number Var; 

CMumber **HHO; // is the vacuum Hamiltonian in the flavour eingenstates 
struct niix{ // contains the mixing matrix 

Number mass [N + 1] ; 

CNumber U[N + 1] [N + 1] ; 
} mixing; 

/* nu are the wave functions, dnu their derivative */ 
CNumber *nu, *dnu; 

int main(int arg, char** argv)-[ 
Number xl=0 . ; 



time_t ti, tf; 
CNumber *onu; 

nu = Cvector (N) ; 

dnu = Cvector (N) ; 
onu = Cvector (N) ; 

ofstream out ("runge.dat"); /* save data in runge.dat */ 

srand (time (&ti)); 
time (&ti) ; 

Var = Varl; 
Number VarO = Var; 

Number h = (VarF - Varl) / INIT_STEPS; 
Number eps = Eps_Error; 
vacuum_ values () ; 



Number x2; 



/* Argument and parameter validation ( 
/* Information output ( 



) */ 
.)*/ 
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/* Main Routines */ 

odeint (nu, N, xl, x2, eps, dist, dist_min, H, deriv, rkqs) ; 
evolute (nu, &out, VarO) ; 

for (int nstp = 1; nstp <= MAX_STEPS; nstp++) 
i /* Take at most MAXSTP steps */ 

for (int i = 1; i <= N; i++) 

onu [i] = nu [i] ; 

vacuuin_values () ; 

if ((VarO + h - VarF) * (VarF - Varl) > 0.0) /* Are we done? */ 
h = VarF - VarO; 
Var = VarO + h; 

odeint (nu, N, xl, x2, eps, dist, 0.0, H, deriv, rkqs); 

if (distance (nu, onu, M) > Prob_Error) 

{ 

//cout « "DECREASE: h=" « h « endl; 

h *= DECREASE; 

Number htemp = ((h < 0) ? 

FMIN (h, -abs ((Varl - VarF) / MAX_STEPS)) : 
FMAX (h, abs ((Varl - VarF) / MAX_STEPS) ) ) ; 
if (htemp != h) 
{ 

evolute (nu, &out, Var); 
VarO = Var; 

} 

h = htemp; 
} 

else 
{ 

evolute (nu, feout, Var); 

h *= INCREASE; 

Number htemp = ((h < 0) ? 

FMAX (h, -abs ((Varl - VarF) / MIN_STEPS)) : 
FMIM (h, abs ((Varl - VarF) / MIN_STEPS) ) ) ; 
if (htemp != h) 
h = htemp; 
VarO = Var; 
} 

if ((Var - VarF) * (VarF - Varl) > 0.0) 

{ /* Are we done? */ 
cout « "t=" « time (&tf) - ti « endl; 
return 0; /* normal exit */ 

} 

} 

cout << "Too many steps in routine evolution_matter ! " << endl; 
return -1; 

} 
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7.2 Example of Hamiltonian definition 

Here is the hamiltonian we use in our sample program: 



/* The Hamiltonian */ 



inline CNumber V(int i, int 

return (i == j == 1) ? 1. : 0.; 

} 

inline CNumber U(int i, int j){ 
return mixing. U[i] [j] ; 

} 

CNumber H (Number r, int i, int j){ 
CNumber HH = 0; 

return HHO[i][j] + V (i, j) * rho (r) * sqrt (2.) * Gf; 

} 



void vacuum_values(){ 
mixing. mass [1] = l.e-2; 
mixing. mass [2] = l.e-1; 
Number thl2 = M_PI / 3. 
Number thl3 = M_PI / 3. 
Number th23 = M_PI / 3. 



nu[l] =1.; 
nu [2] = . ; 

CNumber HO [N + 1] [N + 1] ; 



HHO = matrix ((CNumber) 1 

Number sthl2 = sin (thl2) 

Number cthl2 = cos (thl2) 

Number sthl3 = sin (thl3) 

Number cthl3 = cos (thl3) 

Number sth23 = sin (th23) 

Number cth23 = cos (th23) 



N, 1, N); 



/* 



Three neutrinos 



*/ 



(mixing. U) [1] [1] = CNumber (cthl2 * cthl3) ; 

(mixing. U) [1] [2] = CNumber (sthl2 * cthl3) ; 

(mixing. U) [1] [3] = CNumber (sthl3) ; 

(mixing. U) [2] [1] = CNumber (-sthl2 * cth23 - cthl2 * sth23 * sthl3) ; 
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(mixing. U) [2] [2] = CNumber (cthl2 * cth23 - sthl2 * sth23 * sthl3) ; 
(mixing. U) [2] [3] = CNumber (sth23 * cthl3) ; 

(mixing. U) [3] [1] = CNumber (sthl2 * sth23 - cthl2 * cth23 * sthl3) ; 
(mixing. U) [3] [2] = CNumber (-cthl2 * sth23 - sthl2 * cth23 * sthl3) ; 
(mixing. U) [3] [3] = CNumber (cth23 * cthl3) ; 
break; 
default: 

cerr << "Number of neutrina (" « N « ") not implemented!" « endl; 
exit (-1) ; 

} 

/* 

cout << "Vacuum values :\n"; 
cout << "- HEimiltoniEin: \n" ; 
Hamiltonian HH = H0(); 
cout « HH; 

cout << "- Mixing matrix :\n"; 

cout « * (mixing. U); 

*/ 

for (int i = 1; i <= N; i++) 
for (int j = 1; j <= N; j++) 
{ 

HHO[i] [j] = 0; 

/* HO is the vacuum Hamiltonian in the mass eingenstates */ 
if ((i == 1) && (j == D) 

HO[i] [i] = 1. / pow (10, Var); // This is OK for two neutrina 
else 

HO[i] [j] = 0; 

/* HHO is the vacuum Hamiltonian in the flavour eingenstates */ 

for (int k = 1; k < N; k++) 
for (int 1=1; 1 < N; 1++) 

HHO[i][j] += conj (U (k, i)) * HO[k][l] * U (1, j); 

} 



return; 

} 

The user should provide routines for computing the electron densities at Sun and the Earth 
along the neutrino trajectory. They appear in the matter part of the neutrino hamiltonian. 

7.3 Sample outputs 

The following is the verbatim output of our program 
hamevol.x 

for the values of the parameters which appears in the first information lines (included by 
default in the sample program). 
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./hsmievol.x 

Starting evolution in the Sun 

Used parsmieters: 

MAX_STEPS 100000 1N1T_STEPS 10000 

DECREASE 0.1 INCREASE 5 
Varl -2.39794 VarF -12.3979 

Eps_Error le-08 Prob_Error 0.01 

x2/VarI 0.00881916 x2/VarF 8.81916e+07 



-2. 


,3979 


1 


0.00382 


4. 


,62e- 


-44 


4e-08 


-2. 


,4979 


1 


0.00382 


4. 


,62e- 


-44 


4e-08 


-2. 


,5979 


1 


0.00481 


4. 


,62e- 


-44 


-1.4e-08 


-2. 


,6979 


1 


0.00605 


4. 


,62e- 


-44 


-8.5e-08 


-2. 


,7979 


1 


0.00762 


4. 


,62e- 


-44 


1.5e-09 


-2. 


,8979 


1 


0.00959 


4. 


,62e- 


-44 


-1.7e-08 


-2. 


,9979 


1 


0.0121 


4. 


,62e- 


-44 


3.6e-08 


-3. 


,0979 


1 


0.0152 


4. 


,62e- 


-44 


-2.8e-08 


-3. 


,1979 


1 


0.0191 


4. 


,62e- 


-44 


7.8e-08 
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V CW. IOjUIX^ 




H pcpvi nt.i on 


hamevol-rungekutta . hpp 






MAXSTP 


1000000 


RK algorithm internal parameter 


TINY 


1.0 X 10"^° 


id. 


SAFETY 


0.9 


id. 


PGROW 


-0.2 


id. 


PSHRNK 


-0.25 


id. 


ERRCON 


1.89 X 10"^ 


id. 


hamevol-scunple . cpp 






MAX-STEPS 


1.0 X 10^ 


Steeper method 


MIN-STEPS 


10000 


id. 


INIT-STEPS 


10000 


id. 


DECREASE 


0.1 


id. 


INCREASE 


5.0 


id. 






id. 


hamevol-sample . hpp 






f ermi-MeV 


1.0/197.326 


conversion / — 1/MeV 


m-eV 


fermi-MeV x-^ 


conversion m 1/eV 


Gf 


1.66 X IQ-^^ 


The Fermi constant in 1/eV'^ 


Na 


6.022 X 10^3 


Avogadro number 


RSun 


6.961 X 10^ X m-eV 


Radius of the Sun (1/ eV) 


REarth 


6.378 X lO^x m-eV 


Radius of the Earth (1/eV) 


Eps-Error 


1.0 X 10-8 




Prob-Error 


1.0 X 10-2 




N 


2 


number of equations 


dist 


0.00001 


initial stepsize for Runge-Kutta 


dist-min 


0.0000001 


minimal stepsize for Runge-Kutta 


EARTH 





Program option flag 


SUN 


1 


Program option flag 



Table 2: Here is a list of the most important switches and constants 



Subroutine 


Purpose 


derivs 

runge 

odeint 

rkqs 

rkck 


Computes the derivatives dy/dx 

Given the functions y and their derivatives dy/dx, it returns the advanced solution 
RK driver with adaptive stepsize control. Integrates the starting value over an interval 
with a required accuracy 

Used to monitor accuracy and adjust stepsize during RK integration 

Returns advMnc(Hi solution cm^' mu int.(n'\^al togcMlii^r with the estimate of truucatiou (^Tor 



Table 3: The main subroutines and functions used in the code are reported together with a brief 
explanation of their meaning. 
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